function varargout=fevd_q(varargin)
% This function conducts forcast error variance decomposition
% input
% varargin{1}=gx
% varargin{2}=hx
% varargin{3}=netashock: q shock
% varargin{4}=netashock_all; all shock
% varargin{5}=opt;
%
% output
% rs : fevd result
% 2019/12/15

%--------------------------------------------------------------------------
% input
gx=varargin{1};
hx=varargin{2};
netashock=varargin{3};
netashock_all=varargin{4};
opt=varargin{5};

%--------------------------------------------------------------------------
% choose control variables
gx=gx(opt.var_fevd,:); 

% q shock fevd
nshock=zeros(size(netashock));
nshock(:,1)=netashock(:,1); % use 1st shock

Sig_tmp=zeros(size(gx,1),1); % preallocation
Sig=zeros(size(gx,1),opt.T_fevd);

for h=1:opt.T_fevd % horizon
    
   tmp=gx*hx^(h-1)*nshock;

   % y(t+h)= C1*e(t+h)+C2*e(t+h-1)+...
   % var(y(t+h)|t)
   Sig_tmp=Sig_tmp+diag(tmp*tmp'); 
       
   Sig(:,h)=Sig_tmp; % variable by horizon
      
end

% all shock fevd
nshock=netashock_all;

Sig_tmp=zeros(size(gx,1),1); % preallocation
Sig_all=zeros(size(gx,1),opt.T_fevd);

for h=1:opt.T_fevd % horizon
    
   tmp=gx*hx^(h-1)*nshock;

   % y(t+h)= C1*e(t+h)+C2*e(t+h-1)+...
   % var(y(t+h)|t)
   Sig_tmp=Sig_tmp+diag(tmp*tmp'); 
       
   Sig_all(:,h)=Sig_tmp; % variable by horizon
      
end

% fevd
rs.fevd=Sig./Sig_all;

%--------------------------------------------------------------------------
% output

varargout{1}=rs;




